// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtMetricThresholdFilter.C
// ************************************************************************* //

#include <avtMetricThresholdFilter.h>

#include <vtkCell.h>
#include <vtkCellData.h>
#include <vtkDataSet.h>
#include <vtkDataArray.h>
#include <vtkPointData.h>
#include <vtkPoints.h>
#include <vtkRectilinearGrid.h>
#include <vtkThreshold.h>
#include <vtkUnstructuredGrid.h>

#include <avtCallback.h> 

#include <DebugStream.h>
#include <ImproperUseException.h>
#include <vtkVisItUtility.h>

// ****************************************************************************
//  Method: avtMetricThresholdFilter constructor
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Fri Jun 14 15:42:05 PST 2002
//
// ****************************************************************************

avtMetricThresholdFilter::avtMetricThresholdFilter()
{
}


// ****************************************************************************
//  Method: avtMetricThresholdFilter destructor
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Fri Jun 14 15:42:05 PST 2002
//
//  Modifications:
//
// ****************************************************************************

avtMetricThresholdFilter::~avtMetricThresholdFilter()
{   
}


// ****************************************************************************
//  Method:  avtMetricThresholdFilter::Create
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Fri Jun 14 15:42:05 PST 2002
//
// ****************************************************************************

avtFilter *
avtMetricThresholdFilter::Create()
{
    return new avtMetricThresholdFilter();
}


// ****************************************************************************
//  Method:      avtMetricThresholdFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attsibute object.
//
//  Arguments:
//      a        The attsibutes to use.
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Fri Jun 14 15:42:05 PST 2002
//
// ****************************************************************************

void
avtMetricThresholdFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const MetricThresholdAttributes*)a;
}


// ****************************************************************************
//  Method: avtMetricThresholdFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtMetricThresholdFilter with the given
//      parameters would result in an equivalent avtMetricThresholdFilter.
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Fri Jun 14 15:42:05 PST 2002
//
// ****************************************************************************

bool
avtMetricThresholdFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(MetricThresholdAttributes*)a);
}


// ****************************************************************************
//  Method: avtMetricThresholdFilter::ExecuteData
//
//  Purpose:
//      Sends the specified input and output through the MetricThreshold filter.
//
//  Arguments:
//      in_dr      The input data representation.
//
//  Returns:       The output data representation.
//
//  Programmer: haddox1 -- generated by xml2info
//  Creation:   Fri Jun 14 15:42:05 PST 2002
//
//  Modifications:
//
//    Hank Cihlds, Wed Sep 11 08:50:24 PDT 2002
//    Fixed memory leak.
//
//    Eric Brugger, Wed Jul 30 19:15:53 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

avtDataRepresentation *
avtMetricThresholdFilter::ExecuteData(avtDataRepresentation *in_dr)
{
    //
    // Get the VTK data set.
    //
    vtkDataSet *in_ds = in_dr->GetDataVTK();

    //Get ahold of the data we'll be looking at
    vtkDataArray *scalars = in_ds->GetCellData()->GetScalars();

    if (scalars == NULL)
    {
        scalars = in_ds->GetPointData()->GetScalars();
        if (scalars != NULL)
        {
            avtCallback::IssueWarning("No cell data found, using point data.");

            vtkDataSet *out_ds = ExecuteOnScalarData(in_ds, scalars); 

            avtDataRepresentation *out_dr = new avtDataRepresentation(out_ds,
                in_dr->GetDomain(), in_dr->GetLabel());

            out_ds->Delete();

            return out_dr;
        }
        else
            EXCEPTION0(ImproperUseException);
    }

    vtkUnstructuredGrid *results = vtkUnstructuredGrid::New();

    // We have no idea how many cells we're going to be allocating.
    // But if it's a decent mesh, it shouldn't be huge percentage.
    // In any event, we'll allocate in stages of 256.
    results->Allocate(0, 256);

    // Let's get the points from the input dataset.
    vtkPoints *pts = NULL;
    bool shouldDelete = false;
    switch (in_ds->GetDataObjectType())
    {
    // Easily done, just grab them.
      case VTK_UNSTRUCTURED_GRID:
      case VTK_POLY_DATA:
      case VTK_STRUCTURED_GRID:
        pts = ((vtkPointSet *) in_ds)->GetPoints();
        break;
    
    // If they're a rectilinear grid, we need to get them ourselves
      case VTK_RECTILINEAR_GRID:
        pts = CreateRectilinearPoints((vtkRectilinearGrid *) in_ds);
        shouldDelete = true;
        break;
    
    // If they're any other type of grid, this operator shouldn't be applied
      default:
        EXCEPTION0(ImproperUseException);
    }

    results->SetPoints(pts);
    if (shouldDelete)
        pts->Delete();

    vtkCellData *outcd = results->GetCellData();
    vtkCellData *incd = in_ds->GetCellData();

    outcd->CopyAllocate(incd);

    // We iterate through all the cells in the incoming dataset
    // and run the test
    int i;
    int nCells = in_ds->GetNumberOfCells();
    for (i = 0; i < nCells; i++)
    {
        vtkCell *cell = in_ds->GetCell(i);

        if (PassesTest(cell->GetCellType(), scalars->GetComponent(i,0)))
        {
            // Passes the test, add the cell into our results
            int newid = results->InsertNextCell(cell->GetCellType(),
                                                cell->GetPointIds());
            // And copy over the corresponding cell data.
            outcd->CopyData(incd, i, newid);
        }

    }

    outcd->Squeeze();
    // Move over point data
    results->GetPointData()->PassData(in_ds->GetPointData());    

    avtDataRepresentation *out_dr = new avtDataRepresentation(results,
        in_dr->GetDomain(), in_dr->GetLabel());

    results->Delete();

    return out_dr;
}

// ****************************************************************************
//  Method: avtThresholdFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Indicates the zones no longer correspond to the original problem.
//
//  Programmer: Akira Haddox
//  Creation:   June 20, 2002
//
//  Modifications:
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtMetricThresholdFilter::UpdateDataObjectInfo(void)
{
    GetOutput()->GetInfo().GetValidity().InvalidateZones();
    GetOutput()->GetInfo().GetAttributes().AddFilterMetaData("MetricThreshold");
}

// ****************************************************************************
//  Method: avtThresholdFilter::CreateRectilinearPoints
//
//  Purpose:
//    Find and create a structure for the points of a rectilinear grid.
//
//  Arguments:
//      grid    The rectilinear grid in question.
//
//  Returns:    vtkPoints* structure of the points of the grid. Caller is
//              responsible for deallocating.
//
//  Programmer: Akira Haddox
//  Creation:   June 20, 2002
//
//  Modifications:
//    Kathleen Biagas, Tue Aug 21 16:16:32 MST 2012
//    Preserve coordinate type.
//
// ****************************************************************************
 
vtkPoints *
avtMetricThresholdFilter::CreateRectilinearPoints(vtkRectilinearGrid *grid)
{
    vtkPoints *pts = vtkVisItUtility::NewPoints(grid);

    vtkDataArray *xC = grid->GetXCoordinates();
    vtkDataArray *yC = grid->GetYCoordinates();
    vtkDataArray *zC = grid->GetZCoordinates();

    int numX = xC->GetNumberOfTuples();
    int numY = yC->GetNumberOfTuples();
    int numZ = zC->GetNumberOfTuples();

    pts->SetNumberOfPoints( numX*numY*numZ );

    for (int i = 0; i < numX; i++)
        for (int j = 0 ; j < numY; j++)
            for (int k = 0 ; k < numZ; k++)
            {
                double pt[3];
                pt[0] = xC->GetComponent(i,0);
                pt[1] = yC->GetComponent(j,0);
                pt[2] = zC->GetComponent(k,0);

                pts->SetPoint(k*numX*numY+j*numX+i, pt);
            }
    return pts;
}

// ****************************************************************************
//  Method: avtThresholdFilter::PassesTest
//
//  Purpose:
//      Run the threshold and type test. To pass the test, we must be testing
//      for the type, and the value must not be between the exclusion range.
//
//  Arguments:
//      type       The type of vtkCell.
//      value      The scalar value of the cell.
//
//  Returns:       True on pass (include on output dataset), false on fail.
//
//  Programmer: Akira Haddox
//  Creation:   June 20, 2002
//
//  Modifications:
//    Akira Haddox, Wed Jul  2 08:23:55 PDT 2003
//    Added pixel cell type.
//
// ****************************************************************************

inline bool Between(double _a, double _b, double _c)
{
     return (_a>=_b && _a<=_c);
}

inline bool
avtMetricThresholdFilter::PassesTest(int type, double value)
{
    switch (type)
    {
        case VTK_VOXEL:
        case VTK_HEXAHEDRON:
            return (atts.GetHexahedron() && 
                    !Between(value,atts.GetHex_lower(), atts.GetHex_upper()));
        case VTK_TETRA:
            return (atts.GetTetrahedron() &&
                    !Between(value, atts.GetTet_lower(), atts.GetTet_upper()));
        case VTK_WEDGE:
            return (atts.GetWedge() &&
                    !Between(value, atts.GetWed_lower(), atts.GetWed_upper()));
        case VTK_PYRAMID:
            return (atts.GetPyramid() &&
                    !Between(value, atts.GetPyr_lower(), atts.GetPyr_upper()));
        case VTK_TRIANGLE:
            return (atts.GetTriangle() &&
                    !Between(value, atts.GetTri_lower(), atts.GetTri_upper()));
        case VTK_PIXEL:
        case VTK_QUAD:
            return (atts.GetQuad() && 
                    !Between(value, atts.GetQuad_lower(),
                             atts.GetQuad_upper()));
    }
    return false; //Just in case it's something we haven't considered
}


// ****************************************************************************
//  Method: avtMetricThresholdFilter::ExecuteOnScalarData
//
//  Purpose:
//      Same as normal excecution, but running off of point data instead.
//
//  Arguments:
//      in_ds      The input dataset.
//    scalars       The point scalar data.
//
//  Returns:       The output dataset.
//
//  Programmer: Akira Haddox
//  Creation:   7/2/02
//
//  Modifications:
//    Eric Brugger, Wed Jul 30 19:15:53 PDT 2014
//    Modified the class to work with avtDataRepresentation.
//
// ****************************************************************************

vtkDataSet *
avtMetricThresholdFilter::ExecuteOnScalarData(vtkDataSet *in_ds,
                                              vtkDataArray *scalars)
{
    vtkUnstructuredGrid *results = vtkUnstructuredGrid::New();

    // We have no idea how many cells we're going to be allocating.
    // But if it's a decent mesh, it shouldn't be huge percentage.
    // In any event, we'll allocate in stages of 256.
    results->Allocate(0, 256);

    // Let's get the points from the input dataset.
    vtkPoints *pts = NULL;
    bool shouldDelete = false;
    switch (in_ds->GetDataObjectType())
    {
      // Easily done, just grab them.
      case VTK_UNSTRUCTURED_GRID:
      case VTK_POLY_DATA:
      case VTK_STRUCTURED_GRID:
        pts = ((vtkPointSet *) in_ds)->GetPoints();
        break;
    
      // If they're a rectilinear grid, we need to get them ourselves
      case VTK_RECTILINEAR_GRID:
        pts = CreateRectilinearPoints((vtkRectilinearGrid *) in_ds);
        shouldDelete = true;
        break;
    
      // If they're any other type of grid, this operator shouldn't be applied
      default:
        EXCEPTION0(ImproperUseException);
    }

    results->SetPoints(pts);
    if (shouldDelete)
        pts->Delete();

    vtkCellData *outcd = results->GetCellData();
    vtkCellData *incd = in_ds->GetCellData();

    outcd->CopyAllocate(incd);


    // We iterate through all the cells in the incoming dataset,
    // and run the test
    int i;
    int nCells=in_ds->GetNumberOfCells();
    for ( i = 0; i < nCells; i++)
    {
        vtkCell *cell = in_ds->GetCell(i);

        // Here's where things get a little bit weird
        // What we do is run the test on all the points in the cell, and
        // see if they pass. If any points pass, we'll add the cell into
        // our results. This may be set as an option later, like the Threshold
        // filter, but for now...

        int nPtsCell = cell->GetNumberOfPoints();
        for (int pt = 0; pt < nPtsCell; pt++)
        {
            int globalId = cell->GetPointId(pt);
            if (PassesTest(cell->GetCellType(),
                           scalars->GetComponent(globalId,0)))
            {
                // Passes the test, add the cell into our results
                int newid = results->InsertNextCell(cell->GetCellType(),
                                                    cell->GetPointIds());
                // And copy over the corresponding cell data.
                outcd->CopyData(incd, i, newid);
                break;
            }
        }
    }

    outcd->Squeeze();
    // Move over point data
    results->GetPointData()->PassData(in_ds->GetPointData());    

    return results;
}
